# 사용할 패키지를 추가한다.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
library(dplyr)
library(tidyr)
library(ggmap)
library(tmaptools)
library(raster)
## Warning: package 'raster' was built under R version 4.0.5
library(rgdal)
## Warning: package 'rgdal' was built under R version 4.0.5
library(clValid)
## Warning: package 'clValid' was built under R version 4.0.5
## Warning: package 'cluster' was built under R version 4.0.5
options(scipen=999)

### 버스정류장 clustering

bus <- read.csv("제주특별자치도_버스정류소현황.csv")
bus <-bus[c(3:4)]
head(bus)
##       경도     위도
## 1 126.5609 33.24760
## 2 126.5604 33.24700
## 3 126.5609 33.24751
## 4 126.4631 33.35782
## 5 126.4629 33.35783
## 6 126.4315 33.25145
# # 최적 k 찾기
# 
# b_clvalid <- clValid(bus, 2:20, clMethods="kmeans", validation="internal", maxitems=nrow(bus))
# summary(b_clvalid)


# 위도, 경도 정보
bbox <- c(left=min(bus$경도)-0.05, right=max(bus$경도)+0.05, bottom=min(bus$위도)-0.02, top=max(bus$위도)+0.02)

# 지도 불러오기
map <- get_stamenmap(bbox, zoom=12)

# k=10, 10-means clusterig

b_kmeans <- kmeans(bus,10)

map_bus <- ggmap(map, base_layer=ggplot(bus, aes(x=경도, y=위도), col=b_kmean$cluster))+theme_void()+geom_point(color=b_kmeans$cluster, size=1.5)

map_bus.cluster <- map_bus+geom_point(aes(x=경도, y=위도), data=as.data.frame(b_kmeans$centers),color="blue", fill="blue", size=3)

print(map_bus.cluster)

center <- b_kmeans$centers
center
##        경도     위도
## 1  126.7957 33.50379
## 2  126.5387 33.49376
## 3  126.4349 33.47087
## 4  126.4085 33.26302
## 5  126.8871 33.43497
## 6  126.7684 33.32669
## 7  126.2483 33.28815
## 8  126.2998 33.41706
## 9  126.6616 33.49384
## 10 126.5716 33.26220


### 인구수 map

popul <- read.csv("인구.csv")
popul <-popul[c(4:5)]
popul <- head(popul,10)


pop_bbox <- c(left = 126.08, right=127, bottom=33.1, top=33.6)
p_map <- get_stamenmap(bbox=pop_bbox, zoom=11, maptype="terrain")

map_pop <- ggmap(p_map, base_layer=ggplot(popul, aes(x=경도, y=위도)))+theme_void()+geom_point(color="blue", size=3)

print(map_pop)

center <-rbind(center,popul)
center
##         경도     위도
## 1   126.7957 33.50379
## 2   126.5387 33.49376
## 3   126.4349 33.47087
## 4   126.4085 33.26302
## 5   126.8871 33.43497
## 6   126.7684 33.32669
## 7   126.2483 33.28815
## 8   126.2998 33.41706
## 9   126.6616 33.49384
## 10  126.5716 33.26220
## 11  126.4774 33.48312
## 21  126.5350 33.49716
## 31  126.4969 33.48828
## 41  126.5450 33.47632
## 51  126.3294 33.46201
## 61  126.5383 33.51165
## 71  126.5854 33.52181
## 81  126.6343 33.53438
## 91  126.5655 33.52026
## 101 126.2671 33.41046


### 유동인구수 map

mpopul <- read.csv("유동인구.csv")
mpopul <-mpopul[c(6:7)]
mpopul <- head(mpopul,10)


mpop_bbox <- c(left = 126.08, right=127, bottom=33.1, top=33.6)
mp_map <- get_stamenmap(bbox=mpop_bbox, zoom=11, maptype="terrain")

map_mpop <- ggmap(mp_map, base_layer=ggplot(mpopul, aes(x=경도, y=위도)))+theme_void()+geom_point(color="blue", size=3)

print(map_mpop)

center <-rbind(center,mpopul)
center
##         경도     위도
## 1   126.7957 33.50379
## 2   126.5387 33.49376
## 3   126.4349 33.47087
## 4   126.4085 33.26302
## 5   126.8871 33.43497
## 6   126.7684 33.32669
## 7   126.2483 33.28815
## 8   126.2998 33.41706
## 9   126.6616 33.49384
## 10  126.5716 33.26220
## 11  126.4774 33.48312
## 21  126.5350 33.49716
## 31  126.4969 33.48828
## 41  126.5450 33.47632
## 51  126.3294 33.46201
## 61  126.5383 33.51165
## 71  126.5854 33.52181
## 81  126.6343 33.53438
## 91  126.5655 33.52026
## 101 126.2671 33.41046
## 12  126.5350 33.49716
## 22  126.4774 33.48312
## 32  126.4969 33.48828
## 42  126.3294 33.46201
## 52  126.5450 33.47632
## 62  126.5383 33.51165
## 72  126.5655 33.52026
## 82  126.6343 33.53438
## 92  126.2671 33.41046
## 102 126.3307 33.25741


### 주요장소 최적 k 찾기

place <- read.csv("주요장소.csv")
place <-place[c(1:2)]



# 주요장소 clustering

# 위도, 경도 정보
pl_bbox <- c(left=min(place$경도)-0.05, right=max(place$경도)+0.05, bottom=min(place$위도)-0.02, top=max(place$위도)+0.02)

# 지도 불러오기
pmap <- get_stamenmap(pl_bbox, zoom=12)

# k=10, 10-means clusterig

pl_kmeans <- kmeans(place,10)

map_pl <- ggmap(pmap, base_layer=ggplot(place, aes(x=경도, y=위도), col=pl_kmean$cluster))+theme_void()+geom_point(color=pl_kmeans$cluster, size=1.5)

map_pl.cluster <- map_pl+geom_point(aes(x=경도, y=위도), data=as.data.frame(pl_kmeans$centers),color="blue", fill="blue", size=3)

print(map_pl.cluster)

center <-rbind(center, pl_kmeans$centers)
center
##         경도     위도
## 1   126.7957 33.50379
## 2   126.5387 33.49376
## 3   126.4349 33.47087
## 4   126.4085 33.26302
## 5   126.8871 33.43497
## 6   126.7684 33.32669
## 7   126.2483 33.28815
## 8   126.2998 33.41706
## 9   126.6616 33.49384
## 10  126.5716 33.26220
## 11  126.4774 33.48312
## 21  126.5350 33.49716
## 31  126.4969 33.48828
## 41  126.5450 33.47632
## 51  126.3294 33.46201
## 61  126.5383 33.51165
## 71  126.5854 33.52181
## 81  126.6343 33.53438
## 91  126.5655 33.52026
## 101 126.2671 33.41046
## 12  126.5350 33.49716
## 22  126.4774 33.48312
## 32  126.4969 33.48828
## 42  126.3294 33.46201
## 52  126.5450 33.47632
## 62  126.5383 33.51165
## 72  126.5655 33.52026
## 82  126.6343 33.53438
## 92  126.2671 33.41046
## 102 126.3307 33.25741
## 13  126.4815 33.48367
## 23  126.7408 33.53165
## 33  126.2721 33.27028
## 43  126.5623 33.25835
## 53  126.5474 33.48852
## 63  126.4165 33.47526
## 73  126.8715 33.42594
## 83  126.5192 33.50790
## 93  126.2944 33.42720
## 103 126.5870 33.51688


### 교통량 clustering

car <- read.csv("교통량.csv")
car <-car[c(13:14)]



# 위도, 경도 정보
c_bbox <- c(left=min(car$경도)-0.05, right=max(car$경도)+0.05, bottom=min(car$위도)-0.02, top=max(car$위도)+0.02)

# 지도 불러오기
cmap <- get_stamenmap(c_bbox, zoom=12)

# k=10, 10-means clusterig

c_kmeans <- kmeans(car,10)

map_c <- ggmap(cmap, base_layer=ggplot(car, aes(x=경도, y=위도), col=c_kmean$cluster))+theme_void()+geom_point(color=c_kmeans$cluster, size=1.5)

map_c.cluster <- map_c+geom_point(aes(x=경도, y=위도), data=as.data.frame(c_kmeans$centers),color="blue", fill="blue", size=3)

print(map_c.cluster)

center <-rbind(center, c_kmeans$centers)
center
##         경도     위도
## 1   126.7957 33.50379
## 2   126.5387 33.49376
## 3   126.4349 33.47087
## 4   126.4085 33.26302
## 5   126.8871 33.43497
## 6   126.7684 33.32669
## 7   126.2483 33.28815
## 8   126.2998 33.41706
## 9   126.6616 33.49384
## 10  126.5716 33.26220
## 11  126.4774 33.48312
## 21  126.5350 33.49716
## 31  126.4969 33.48828
## 41  126.5450 33.47632
## 51  126.3294 33.46201
## 61  126.5383 33.51165
## 71  126.5854 33.52181
## 81  126.6343 33.53438
## 91  126.5655 33.52026
## 101 126.2671 33.41046
## 12  126.5350 33.49716
## 22  126.4774 33.48312
## 32  126.4969 33.48828
## 42  126.3294 33.46201
## 52  126.5450 33.47632
## 62  126.5383 33.51165
## 72  126.5655 33.52026
## 82  126.6343 33.53438
## 92  126.2671 33.41046
## 102 126.3307 33.25741
## 13  126.4815 33.48367
## 23  126.7408 33.53165
## 33  126.2721 33.27028
## 43  126.5623 33.25835
## 53  126.5474 33.48852
## 63  126.4165 33.47526
## 73  126.8715 33.42594
## 83  126.5192 33.50790
## 93  126.2944 33.42720
## 103 126.5870 33.51688
## 14  126.3975 33.28989
## 24  126.2826 33.25185
## 34  126.2001 33.34196
## 44  126.6102 33.28166
## 54  126.8047 33.42854
## 64  126.2464 33.28657
## 74  126.3675 33.40089
## 84  126.6056 33.48236
## 94  126.4593 33.45413
## 104 126.2829 33.39465


options(scipen=999)

### 각 센터들을 다시 clustering

center
##         경도     위도
## 1   126.7957 33.50379
## 2   126.5387 33.49376
## 3   126.4349 33.47087
## 4   126.4085 33.26302
## 5   126.8871 33.43497
## 6   126.7684 33.32669
## 7   126.2483 33.28815
## 8   126.2998 33.41706
## 9   126.6616 33.49384
## 10  126.5716 33.26220
## 11  126.4774 33.48312
## 21  126.5350 33.49716
## 31  126.4969 33.48828
## 41  126.5450 33.47632
## 51  126.3294 33.46201
## 61  126.5383 33.51165
## 71  126.5854 33.52181
## 81  126.6343 33.53438
## 91  126.5655 33.52026
## 101 126.2671 33.41046
## 12  126.5350 33.49716
## 22  126.4774 33.48312
## 32  126.4969 33.48828
## 42  126.3294 33.46201
## 52  126.5450 33.47632
## 62  126.5383 33.51165
## 72  126.5655 33.52026
## 82  126.6343 33.53438
## 92  126.2671 33.41046
## 102 126.3307 33.25741
## 13  126.4815 33.48367
## 23  126.7408 33.53165
## 33  126.2721 33.27028
## 43  126.5623 33.25835
## 53  126.5474 33.48852
## 63  126.4165 33.47526
## 73  126.8715 33.42594
## 83  126.5192 33.50790
## 93  126.2944 33.42720
## 103 126.5870 33.51688
## 14  126.3975 33.28989
## 24  126.2826 33.25185
## 34  126.2001 33.34196
## 44  126.6102 33.28166
## 54  126.8047 33.42854
## 64  126.2464 33.28657
## 74  126.3675 33.40089
## 84  126.6056 33.48236
## 94  126.4593 33.45413
## 104 126.2829 33.39465
# 최적 k 찾기

ct_clvalid <- clValid(center, 2:8, clMethods="kmeans", validation="internal", maxitems=nrow(center))
summary(ct_clvalid)
## 
## Clustering Methods:
##  kmeans 
## 
## Cluster sizes:
##  2 3 4 5 6 7 8 
## 
## Validation Measures:
##                            2       3       4       5       6       7       8
##                                                                             
## kmeans Connectivity   6.6524  9.4627 11.5246 20.2984 23.3496 28.7849 34.3456
##        Dunn           0.0418  0.3076  0.3569  0.0831  0.0831  0.0831  0.1282
##        Silhouette     0.5206  0.5543  0.5986  0.5776  0.5481  0.5594  0.5103
## 
## Optimal Scores:
## 
##              Score  Method Clusters
## Connectivity 6.6524 kmeans 2       
## Dunn         0.3569 kmeans 4       
## Silhouette   0.5986 kmeans 4
# 위도, 경도 정

ct_bbox <- c(left = 126.08, right=127, bottom=33.1, top=33.6)
ctmap <- get_stamenmap(bbox=ct_bbox, zoom=11, maptype="terrain")


# k=4, 4-means clusterig

ct_kmeans <- kmeans(center,4)

map_ct <- ggmap(ctmap, base_layer=ggplot(center, aes(x=경도, y=위도), col=ct_kmean$cluster))+theme_void()+geom_point(color=ct_kmeans$cluster, size=1.5)


map_ct.cluster <- map_ct+geom_point(aes(x=경도, y=위도), data=as.data.frame(ct_kmeans$centers),color="blue", fill="blue", size=3)

print(map_ct.cluster)

ct_kmeans$centers
##       경도     위도
## 1 126.2975 33.33642
## 2 126.5695 33.50498
## 3 126.4400 33.47507
## 4 126.7347 33.38376